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Abstract 

Any model of tides is based on a specific hypothesis of how the geometric lag <5i depends 
on the tidal-flexure frequency X- Some authors put 5\ ~ x°> the others set <5i ~ x ■ The 
actual dependence determined by the planet's rheology. A particular form of this depen- 
dence will fix the form of the frequency dependence of the tidal quality factor Q . Since at 
present we know the shape of the function we can reverse our line of reasoning and 

single out the appropriate frequency-dependence of the lag. The dependence is different 
■ from those employed hitherto, and makes a considerable alteration in the theory of tides. 

The second alteration is needed to evade some difficulties inherent in Kaula's model. 



Kaula's expansion is divergent for eccentricities higher than a certain value. It is also 
subject to "Goldreich's admonition": due to nonlinearity, the quality factors introduced 



for higher harmonics are badly defined (Goldreich 1963). To avoid these predicaments, we 
abandon the expansion and introduce an instantaneous tidal frequency that bears a de- 
pendence upon the moon's true anomaly (concept prompted by the WKB approximation 
in quantum mechanics). We also introduce an overall quality factor, the only Q factor 
emerging in the model. This factor is a function of the instantaneous tidal frequency and 
therefore it, too, depends on the true anomaly. 

The third alteration will be our departure from the popular belief that Q is inversely 
proportional to the "tangential lag" (the angle 8\ subtended at the planet's centre between 
the moon and the bulge). That this is correct only in the limit of a circular equatorial 
orbit was pointed out by Kaula (1968) - see the caveat preceding his derivation of formula 
(4.5.19) for the damping rate. Generalisation of Kaula's formula to a finite eccentricity 
and inclination yields extra terms. These become leading on crossing the synchronous 
orbit, when the "tangential lag" vanishes. We show, by some different method, that in 
the general case Q is interconnected, in the standard way, not with the "tangential lag" , 

I 2 2~ 

but with the overall lag 5 = y (tangential lag) + (radial lag) . While the "tangential 
lag" 6i originates because the location of the bulge lags, the "radial" lag comes into play 
since the height of the bulge lags too. 

It turns out that, for an arbitrary eccentricity and inclination, the total lag and the 
instantaneous x an d Q are always interconnected via the same relations as in the trivial 
case of an equatorial circular orbit. This makes our model concise and self-explanatory. 

The model addresses only the land tides, and therefore is intended to be a tool for 
exploring the dynamics of the Martian satellites. It is applicable to the Earth-Moon 
system only for the aeons preceding the formation of oceans. 
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1 The Gerstenkorn-MacDonald-Kaula and 
Singer-Mignard theories of bodily tides 

If a satellite is located at a planetocentric position r, it generates a tidal bulge that either 
advances or retards the satellite motion, depending on the interrelation between the planetary 
spin rate u p and the tangential part of satellite's velocity v divided by r = It is convenient 
to imagine that the bulge emerges beneath a fictitious satellite located at 

r f = r+f , (1) 

— *■ 

where the position lag / is given by 

/ = At (£ p xr - v ) , (2) 

At being the time lag between the real and fictitious tide-generating satellites. We shall also 
introduce the angular lag as 

r I /I A* _ 

o = 1 — - = — \ujp x r — v\ . (3) 



For zero eccentricity, S is simply the absolute value of the angle Si subtended at the planet's 
centre between the satellite and the tidal bulge, as on Fig. [IJ For non-circular orbits, our S 
differs from the absolute value of the subtended angle Si . We have deliberately arranged for 
this difference, to make our S embrace both the "tangential" and the "radial" lagging, i.e., 
the lag in the position and also that in the height of the bulged] 

The imaginary satellite is merely a way of illustrating the time lag between the tide-raising 
potential and the distortion of the body. This concept implies no new physics, and is but 
a convenient figure of speech employed to express the fact that at each instance of time the 
dynamical tide is modelled with a static tide where all the time-dependent variables are shifted 
back by At, i.e., (a) the moon is rotated back by v At , and (b) the attitude of the planet 
is rotated back by uj p At . From the viewpoint of a planet-based observer, this means that a 
dynamical response to a satellite located at r is modelled with a static response to a satellite 
located at f/ = r — At (v — uj p x r) . 

This delayed response implies a lag in the location of the bulge, and a lag in its height. As 
a result, the tidal dissipation consists of two inputs - the "tangential" and the "radial" one. 

Our S defined by ([3]) is equal to y (tangential lag) + (radial lag) . This, together with the 
interconnection between our S and the overall quality factor, ensures that our model indeed 
includes both lags - that in the position and that in the height of the tide. Leaving details 
for Section 3 below, we would emphasise already at this point that the well-known relation 
l/Q = tan 25 interconnects the quality factor not with the "tangential lag" Si but exactly 
with the total lag that includes both the "tangential" and "radial" part. 

The earliest efforts aimed at modeling bodily tides were undertaken in the end of the XIX th 
century by Darwin (1879, 1880). Furthering of this line of research by Gerstenkorn (1955), 

The necessity of including the "radial" part has long been pointed out in the literature - for example, in 
the book by Murray and Dermott (1999), pp. 170 - 171, this issue is discussed in regard to dissipation inside 
satellites. We consider this effect in regard to the dissipation within the planet. 
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Figure 1: The position lag / and the angular lag 5\ of a satellite located below the synchronous 
orbit. (For a moon located above the synchronous orbit both lags will be pointing in an opposite 
direction.) Be mindful that 5 defined through ([3]) does not coincide with the subtended angle 
5\ , except for circular orbits. 
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MacDonald (1964) and Kaula (1964) gave birth to what we now may call 'the old theory of 
tides'. In their works, Gerstenkorn, MacDonald, and KaukU assumed the subtended angle to 
be a fixed constant. With the latter assumption regarded as a critical deficiency, their theory 
soon fell into disuse in favour of the approach offered by Singer and Mignard. Singer (1968) 
suggested that the subtended angle should be proportional to the principal frequency x of the 
tide. This was equivalent to setting At constant in (j2]-[2]). Singer applied this new theory 
to the Moon and to Phobos and Deimos. Sadly, Singer's interesting work is seldom cited. A 
detailed mathematical development of Singer's idea can be found in Mignard (1979, 1980) who, 
actually, completely avoided using the lag angle and operated only with the position and time 
lags. Later, Singer's assumption of a constant At was employed also by Touma & Wisdom 
(1994) and Peale & Lee (2000). 



2 The tidal frequency \ an d the quality factor Q 
2.1 How to define the tidal frequency? 

Singer (1968) approximated the principal frequency of tidal flexure, Xi with the expression 
2 \n — uj p \ , where n and uo p = | uj p \ stand for the satellite's mean motion and the planet's spin 
rate. Thus he restricted the applicability realm of his theory to the case of low inclinations 
and eccentricities. While it is clear that, for an equatorial circular orbit, the satellite velocity 
relative to the surface is 



the tidal frequency is 



and the angular lag is 



oj p x r — v | = r | uj p — n | , (4) 



X = 2 | Up - n | , (5) 



6 = — \UpXr - v\ = — x , (6) 



the question remains how can this machinery be extended to the generic case when neither 
% nor e is small. In particular, how can one define the tidal frequency (or, perhaps, the 
principal tidal frequency and the higher frequencies) when the satellite may appear over a 
different point on the surface after each revolution, thus making the entire notion of flexure 
cycle hard to define? 



2.2 Goldreich's admonition: 

a general difficulty stemming from nonlinearity 

Generalisation of the above construction was offered already by Darwin (1908) and, later, by 
Jeffreys (1961) and Kaula (1964). These works' starting point was that each elementary vol- 
ume of the planet is subject to a tide-raising potential, which in general is not periodic but can 
be expanded into a sum of periodic terms. They then employed the linear approximation intro- 
duced by Love, according to which the tidal perturbations of the potential yield linear response 

2 To be more exact, Kaula's starting point was an assertion of constancy of the tidal quality factor. 
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of the shape and linear variations of the stress. In extension of the linearity approximation, 
these authors also implied that the overall dissipation inside the planet may be represented as 
a sum of attenuation rates corresponding to each periodic disturbance: 



( . . . ) standing for averaging over flexure cycle, E(x i ) denoting the energy of deformation 
at the frequency Xi > an d Q{Xi) being the quality factor of the material at this frequency. 
Introduced empirically as a means to figleaf our lack of knowledge of the attenuation process 
in its full complexity, the notion of Q has proven to be practical due to its smooth and universal 
dependence upon the frequency and temperature. At the same time, this empirical treatment 
has its predicaments and limitations. Its major inborn defect was brought to light by Peter 
Goldreich who pointed out that the attenuation rate at a particular frequency depends not only 
upon the appropriate Fourier component of the stress, but also upon the overall stress. This 
happens because for real minerals each quality factor bears dependence not only on the 

frequency x% > but also on the Xi component of the stress and, most importantly, also on the 
overall stress. This, often-neglected, manifestation of nonlinearity may be tolerated only when 
the amplitudes of different harmonics of stress are comparable. However, when the amplitude 
of the principal mode is orders of magnitude higher than that of the harmonics (tides being 
the case), then the principal mode will, through this nonlinearity, make questionable our entire 
ability to decompose the overall attenuation into a sum over frequencies. Stated differently, 
the quality factors corresponding to the weak harmonics will no longer be well defined physical 
parameters. 

Here follows a quotation from Goldreich (1963): 
"... Darwin and Jeffreys both wrote the tide-raising potential as the sum of periodic potentials. 
They then proceeded to consider the response of the planet to each of the potentials separately. 
At first glance this might seem proper since the tidal strains are very small and should add 
linearly. The stumbling block in this procedure, however, is the amplitude dependence of the 
specific dissipation function. In the case of the Earth, it has been shown by direct measurement 
that Q varies by an order of magnitude if we compare the tide of frequency 2u — 2n with 
the tides of frequencies 2u — n , 2u — 3n , and . This is because these latter tides have 
amplitudes which are smaller than the principle tide ( of frequency 2uj — 2n ) by a factor of 
eccentricity or about 0. 05. It may still appear that we can allow for this amplitude dependence 
of Q merely by adopting an amplitude dependence for the phase lags of the different tides. 
Unfortunately, this is really not sufficient since a tide of small amplitude will have a phase lag 
which increases when its peak is reinforcing the peak of the tide of the major amplitude. This 
non-linear behaviour cannot be treated in detail since very little is known about the response of 
the planets to tidal forces, except for the Earth. " 

On these grounds, Goldreich concluded the paragraph with an important warning that we 
"use the language of linear tidal theory, but we must keep in mind that our numbers are really 
only parametric fits to a non-linear problem. " 

In order to mark the line beyond which this caveat cannot be ignored, let us first of all 
recall that the linear approximation remains applicable insofar as the strains do not approach 



(7) 



where, at each frequency Xi 
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the nonlinearity threshold, which for most minerals is of order 10~ 6 . On approach to that 
threshold, the quality factors may become dependent upon the strain magnitude. In other 
words, in an attempt to extend the expansion ([7| - |Sj) to the nonlinear case, we shall have 
to introduce, instead of Q(Xi) > some new functions Q{Xi, E peak (Xi) , E overaU ) . (Another 
complication is that in the nonlinear regime new frequencies will be generated, but we shall 
not go there.) Now consider a superposition of two forcing stresses - one at the frequency X\ 
and another at \ 2 . Let the amplitude E peak (\\ ) be close or above the nonlinearity threshold, 
and E Bak {x 2 ) be by an order or two of magnitude smaller than E eak (Xi) ■ To adapt the 
linear machinery ([7] - [8]) to the nonlinear situation, we have to write it as 

E peak \Xl) E peak (X 2 ) 



{E) {El) + {E2) Xl Q(xJ a E peak (Xi)) X " Q(x 2 ,E peak ( Xx ), E pe Jx 2 )) ' 

the second quality factor bearing a dependence not only upon the frequency X2 and the 
appropriate magnitude E eak (x 2 ) > but also upon the magnitude of the first mode, E eofc (Xi) > 
- this happens because it is the first mode which makes a leading contribution into the overall 
stress. Even if ([9]) can be validated as an extension of (JT1 - E|) to nonlinear regimes, we should 
remember that the second term in ([9]) is much smaller than the first one (because we agreed 
that E peak (x 2 ) ^ E peak {x\) )• This results in two quandaries. The first one (not mentioned 
by Goldreich) is that a nonlinearity-caused non-smooth behaviour of Q{X\ j E peak (Xi) ) wn l 
cause variations of the first term in (Q, which may exceed or be comparable to the entire 
second term. The second one (mentioned in the afore quoted passage from Goldreich) is the 
phenomenon of nonlinear superposition, i.e., the fact that the smaller-amplitude tidal harmonic 
has a higher dissipation rate (and, therefore, a larger phase lag) whenever the peak of this 
harmonic is reinforcing the peak of the principal mode. Under all these circumstances, fitting 
experimental data to ([9]) will become a risky business. Specifically, it will become impossible 
to reliably measure the frequency dependence of the second quality factor; therefore the entire 
notion of the quality factor will, in regard to the second frequency, become badly defined. 

The admonition by Goldreich had been ignored until an alarm sounded. This happened 
when the JPL Lunar-ranging team applied the linear approach to determining the frequency- 
dependence of the Lunar quality factor, Q(x) ■ They obtained (or, as Peter Goldreich rightly 
said, fitted their data to) the dependency Q(x) ~ X a (Williams et al. 2001). The value of 
the exponential for the Moon turned out to be negative: a = — 0.07 , a result firmly tabooed 
by the condensed-matter physics for this range of frequencies (Karato 2007). 

One possible approach to explaining this result may be the following. Let us begin with a 
very crude estimate for the tidal strains. For the Moon, the tidal displacements are of order 
0.1m (zero to peak). The most rough estimate for the strain can be obtained through dividing 
the displacement by the radius of the body. While this ratio is almost twenty times less than 
the 10 -6 nonlinearity threshold, we should keep in mind that in reality the distribution of the 
tidal strain is a steep function of the radius, with the strain getting its maximum near the 
centre of the body, as can be seen from equations (48.17) in Sokolnikoff (1956). The values 
of strain can vary in magnitude, over the radius, by about a factor of fivejfl Another order 
of magnitude, at least, will come from the fact that the Moon is inhomogeneous and that its 
warmer layers are far more elastic than its rigid surface. As a result, deep in the Lunar interior 
the tidal strain will exceed the afore mentioned nonlinearity threshold. Hence the uncertainties 
in determination of Q(x) by the Lunar-ranging team. 



3 Accordingly, the tidal-energy radial distribution, too, is strongly inhomogeneous, with a maximum in the 
planet's centre - see Fig. 3 in Peale & Cassen (1978). 
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We shall not dwell on this topic in quantitative detail, leaving it for a future work. Our 
only goal here has been to draw the readers' attention to the existing difficulty stemming from 
the shortcomings of the extension of ([7] - [8]) to nonlinear regimes. 

2.3 Specific problems with the Kaula expansion 

The effect of lags in the tidal bulge can be described by expanding the potential of the tidally 
distorted planet in zonal harmonics, each of which is then transformed into a sum of spherical 
ones, with the further intention of expressing each term via the orbital elements and mean 
anomalies of the tide-raising and tide-disturbed satellites. A comprehensive development in this 
direction carried out by Kaula (1964) was later employed and quoted in numerous papers and 
books - see, for example, section VI in Lambeck (1980). However, neither Kaula nor Lambeck 
adumbrated the applicability realm of this approach, nor did they dwell on its drawbacks. 
Below we shall mention several limitations, which restrict the usability of this method. 

First of all, this theory is subject to Goldreich's admonition. Indeed, in the Kaula ex- 
pansion each frequency has its own phase lag, which must be interconnected with the quality 
factor appropriate to this frequency. For the reason described in the previous subsection, this 
inconvenience may be circumvented only in the linear case. In the nonlinear situations (like, 
say, dissipation inside Phobos and, most likely, inside the Moon), this becomes a problem, 
because the quality factors of higher harmonics become badly defined when the magnitudes of 
these harmonics are much lower than that of the principal mode0 

The second problem with Kaula's series is that the derivation of the spherical-harmonics 
expansion - formula (1) in Kaula (1964) or formula (6.3.1) in Lambeck (1980) - contains 
a transition from the true anomaly v to the mean anomaly M , transition based on the 
trigonometric expansion of v — M as a Fourier sine series in M . This expansion turns out to 
be convergent only for e < 0.6627434 . (Murray and Dermott 2000) As a result, the entire 
spherical-harmonics machinery remains legitimate only below this value of e . One negative 
consequence of this mishap is that the theory of Kaula (1964) cannot be applied to check 
the possibility of tidal capture hypothesised by Singer (1965) - exploration of this mechanism 
makes it necessary to study tidal effects at initial eccentricities close to unity. 

2.4 Quest for a new approach 

In our opinion, the right way out of the afore described predicament lies not along the paths of 
mathematical rigour but through a careful choice of the particular physical meaning wherewith 
we want to endow the term "tidal frequency." 

In the theory of tides, the notion of principal frequency is employed to make use of the 
attenuation theory borrowed from acoustics. There the basic fact is that, for a sample of 
material (an important clause to reject later), the time-averaged rate of energy losses produced 
by small-amplitude alternating stresses can be expanded over the frequencies Xi involved: 



4 Goldreich (1966) employed the Kaula expansion to explore the Lunar-orbit history. To circumvent his 
own admonition, Goldreich had to assume all the lags to be the same, and to be constant. This, however, was 
long before the actual low-frequency behaviour of Q was established. 




(10) 
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where 

E (y.) 

E(y.) = - v. v™ kKA ' lJ . (11) 

As the average dissipation rate (E) is, by definition, the energy damped over a cycle, divided 
by the period 2ir/x i , then one can shape the expression for the quality factor into the form: 



AE , (v.) = -2tt (12) 



The advantage of definition (fT2l) is that it operates with only one cycle of flexure. (Even a half- 
cycle would be sufficient.) In no way does this definition refer to the afore emphasised clause 
about a certain sample of material. Rather, (fT2l) may be applied to a wave train and even 
to a solitary wave travelling through a medium; in this situation no fragment of the medium 
experiences a repeated pattern of strain. The solitary wave travels through this fragment just 
once. Nevertheless, ( fT2|) still enables us to talk about the attenuation rate in regard to this 
wave. Let us apply this reasoning to tides and, specifically, to the situation when the satellite 
appears, after each revolution, over a different point of the surface. It is now unimportant that 
the stress in a particular fragment of the Earth material lacks periodicity, and that we cannot 
clearly define a flexure frequency for a particular fragment. As agreed, we shall apply ffT2]) not 
to a fixed fragment, but to the travelling tidal bulge considered as a solitary half- wave of a 
frequency equal to the velocity of the satellite relative to the ground beneath it, divided by r : 

a = - \uj p x r — v\ , (13) 
r 

Had the satellite's pull created only one elevation on the planet's surface, we would employ ([HI) 
to write E = -aE peak /Q(a) or, equivalently, AE cyc[e (a) = -2nE peak /Q(a) . However, since 
the moon's gravity produces two uplifts - one on the facing side and another on the opposite 
side of the planet - then the actual dissipated energy will be doubled: E = —2crE peak /Q(a) 
or, the same, AE e (2o) = —2irE peak /Q(cr) . This can be shaped into the form similar to fj 1 XI) 



or into the form like f[T2|) : 



cycle 

where 



AE , (x) = - 2 7T ^4 (15) 

Q(x) 



X = 2 a = - I uL x r — v I , (16) 
r 

is effectively playing the role of the tidal frequency. This is its definition: in our model, the 
tidal frequency is the quantity x > insertion whereof in (|H|) or (|T2|) gives the correct dissipation 
rate. It remains to be said that in f[T5l-[T6l) we slightly cheated by substituting Q(a) = Q(x/2 
with Q(x) ■ Formally, this trick can be justified as a convenient redefinition of the Q factor 



5 In reality, this formal justification is, of course, redundant, because the dependence Q(x) is very slow, 
and noticeable changes of Q take place over frequencies ranges that span, at least, an order of magnitude. In 
realistic seismological measurements, an increase of the frequency by some factor of two yields a variation in 
the value of Q , that will barely make it over the error bar. 
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The so-defined frequency can be expressed via the planet's spin rate u p , the true anomaly v 
of the satellite, and the orbital elements. Denoting the radial component of the moon's velocity 
with v r , and its planetocentric distance with r , we write (see the Appendix for details): 

| lj p x r — v\ = r \J {u — Up cosi) 2 + u 2 sin 2 i cos 2 (u + v) + v 2 r~ 2 , (17) 

whence it can be derived (once again, see the Appendix) that 

X = 2 \ l{y — ujp cosz) 2 + uj 2 sin 2 i cos 2 (cj + v) + v 2 r~ 2 



\ 



n 



(1 + e cos z/) z 



1 - e 2 ) 



2^/2 



Up COS I 



o ■ 2 • 9/ \ (1 + e cosz/) . 2 

u z sin i cos^fu; + v) + n z e z ■ "" 



;i - e 2 y 



Sill V 



Here n , e , % , u stand for the moon's mean motion, eccentricity, inclination, and the argu- 
ment of the pericentre, while u p is the rotation rate of the planet. 



Let us draw conclusions: 



1. The simple approximation (jU - [6]) can be used only at a relatively late stage of the 
satellite's life - when its orbit is well circularised and damped down toward the equator. 

2^ In the generic case of inclined and eccentric orbits, (jU - [5]) must be substituted with 
ffTTl- [T8l) . while the expression ([6]) for the angular lag will preserve its form: 

5 = At - \C3pXr - v\ = ^ X , (19) 
r 2 

X being now given by (fT8l) . The generalised formalism (fl8l - [T9l) will render the satellite orbit 
evolution right after the capture, when the eccentricity (and, possibly, the inclination) is high. 

3. In our model, the above definition for the angular lag 5 has been deliberately chosen 
to differ, for non-circular orbits, from the angle 5\ subtended at the planet's centre between 
the satellite and the tidal bulge, as on Fig. 1. Had we not arranged for this difference, our 
model would be taking into account only the "tangential" part of the tidal dissipation (one 
associated with the tangential motion of the satellite) § Our definition (TT9T) ensures that both 
the "tangential" and "radial" parts of the tidal dissipation contribute to the lag 5 and, thereby, 
to the overall quality factor Q interconnected with S through ( !20l) . So our model contains a 
radial- velocity-squared term in ([171 - fT8l) . a term that prevents X from passing through zero 
at resonance crossings. 



6 That such a truncated model would be not only incomplete but also violently unphysical can be seen from 
the fact that, on each passage of a satellite through a spin-orbit resonance, the truncated model would cause 
the torque to formally approach infinity and then to suddenly reverse to the infinity of an opposite sign. Indeed, 
since S is proportional to x to a negative power, then in the absence of the "radial" dissipation \ would be 
passing through zero on each crossing of the synchronous orbit. To spare the model from such vicissitudes, we 
must define the lag as in ([TO)) . 
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4. Our model departs from the traditional approach on one more crucial point. Taking 
into consideration the Goldreich admonition, we refrain from decomposing the overall tide into 
modes and from attributing a separate quality factor to each of these. Instead, we introduce 
only one tidal frequency that bears an explicit dependence upon the satellite's true anomaly, 
and corresponds to an appropriate overall quality factor. 

The concept of an instantaneous frequency as a function of the true anomaly is an evident 
analogue to the evolving phase in the WKB approximation in quantum mechanics. This 
mathematical analogy, though, should not be taken too literally, because in quantum mechanics 
this method is used to solve a linear Schrodinger equation, while here we employ it in a highly 
nonlinear and dissipative context. 



3 The quality factor Q and the geometric lag angle 5. 



— Mr. Wallace, are you OK? 

— No, I am very far from being OK... 

"Pulp Fiction" 

During tidal flexure, the energy attenuation through friction is, as ever, accompanied by 
a phase shift between the action and the response. The tidal quality factor is interconnected 
with the phase lag e and the angular lag 5 via 

Q' 1 = tane = tan 26 (20) 

or, for small lag angles, 

Q- 1 « e = 2 5 . (21) 

The doubling of the lag is a nontrivial issue. Many authors erroneously state that Q~ l is equal 
simply to the tangent of the lag, with the factor of two omittedJll In Appendix 2.1, we offer a 
simple illustrative calculation, which explains whence this factor of two stems. 

As an aside, we would mention that Zschau (1978) contested the entire possibility of ex- 
pressing the quality factor via only one lag. His point stemmed from the fact that attenuation 
of land tides, in distinction from earthquakes, takes place both due to rigidity and self-gravity 
of the body. In Appendix A. 3 below we explain why a typical terrestrial planet or small body 
will be exempt from Zschau's objection, so that formulae (1201 - 121]) will hold for them. 

7 Rainey & Aharonson (2006) assume that Q _1 is equal to the tangent of the geometric lag. As a result, they 
arrive to a value of Q that is about twice larger than those obtained by the other teams, for example by Bills 
et al (2005). As pointed out by Valery Lainey, in the work by Bills et al (2005) one letter, 7 , is used to denote 
two different angles. In their paper, prior to equation (24), 7 signifies the geometric lag (in our notations, 
61 ). Further, in their equations (24) and (25), Bills et al employ the notation 7 to denote the phase lag (in 
our notations, e, which happens to be equal to 2 6%). With this crucial caveat, Bills' equation Q = l/tan7 
is correct. This mess in notations has not prevented Bills et al (2005) from arriving to a reasonable value of 
the Martian quality factor, 85.58 ± 0.37 . (A more recent study by Lainey et al (2006) has given a comparable 
value of 79.91 ± 0.69 .) 
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The popular oversight interconnects the quality factor with the "tangential lag," i.e., with 
the geometric angle 8\ depicted on Fig. 1. In our terms, this would read as Q^ 1 = tan25i or 
Q^ 1 « 2 81 . While correct in the limit of a circular equatorial orbit (see Appendix 2.1 below), 
in the general case this relation is plainly wrong. Nonetheless it was employed by many. Kaula 
(1968) explicitly acknowledged that his proof of this relation is valid only in the said limit^l To 
treat the case of an elongated inclined orbit, one can amend Kaula's formula with the terms 
he dropped. In Appendices A. 2. 2 - A. 2. 3 below we shall offer an easier development, from 
which it will be seen that, for an eccentric inclined orbit, the "tangential lag" 5i should be 
substituted with the total lag 5 = \J (tangential lag) 2 + (radial lag) 2 defined through ([3]). 
This will explain why in our formulae ([20] - EJJ) we use the total lag 5 , not the subtended angle 
Si. 

Taken the above elucidation, formulae (12131 - l2~Tj) look reasonable: the higher the quality 
factor, the lower the damping rate and, accordingly, the smaller the lag. What look very 
far from being OK are the frequency dependencies ensuing from the assertions of 5 being 
either constant or linear in frequency: the Gerstenkorn-MacDonald-Kaula theory implies that 
Q ~ x° > while the Singer-Mignard theory yields Q ~ x 1 ? neither option being in agreement 
with the geophysical data. 

4 Dissipation in the mantle. 

4.1 Generalities 

Back in the 60s and 70s of the past century, when the science of low-frequency seismological 
measurements was yet under development, it was widely thought that at long time scales the 
quality factor of the mantle is proportional to the inverse frequency. This fallacy proliferated 
into planetary astronomy where it was received most warmly, because the law Q ~ 1/x 
turned out to be the only model for which the linear decomposition of the tide gives a set of 
bulges displaced from the direction to the satellite by the same angle. Any other frequency 
dependence Q(x) entails superposition of bulges corresponding to the separate frequencies, 
each bulge being displaced by its own angle. This is the reason why the scaling law Q ~ 1/x , 
long disproved and abandoned in geophysics (at least, for the frequency band of our concern), 
still remains a pet model in celestial mechanics of the Solar system. 

Over the past twenty years, a considerable progress has been achieved in the low-frequency 
seismological measurements, both in the lab and in the field. Due to an impressive collective 
effort undertaken by several teams, it is now a firmly established fact that for frequencies down 
to about ~ 1 yr^ 1 the quality factor of the mantle is proportional to the frequency to the power 
of a positive fraction a . This dependence holds for all rocks within a remarkably broad band 
of frequencies: from several MHz down to about 1 yr -1 . 

At timescales longer than 1 yr, all the way to the Maxwell time (about 100 yr), attenuation 
in the mantle is defined by viscosity, so that the quality factor is, for all minerals, well approx- 
imated with 7]x/M , where rj and M are the sheer viscosity and the sheer elastic modulus of 

8 On page 201 of his book, Kaula (1968) calculates the dependence of E upon the tangential lag (which he 
calls (5 , and which we term Si - mind the difference in notations). Kaula rightly prefaces his derivation with 
the words "If the eccentricity and inclination of the moon's orbit are neglected", because his formula (4.5.19) 
is correct in this approximation only. Generalisation of Kaula's formula yields extra terms dependent upon the 
inclination and eccentricity. These terms become leading on approach to or crossing of the synchronous orbit. 
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the mineral. Although the values of both the viscosity coefficients and elastic moduli greatly 
vary for different minerals and are sensitive to the temperature, the overall quality factor of 
the mantle still scales linear in frequency. 

There still is no consensus in the seismological community in regard to the time scales 
exceeding the Maxwell time. One viewpoint (incompatible with the Maxwell model) is that 
the linear law Q ~ x extends all the way down to the zero- frequency limit (Karato 2007). 
An alternative point of view (prompted by the Maxwell model) is that at scales longer than 
the Maxwell time we return to the inverse-frequency law Q ~ l/x- 

All in all, we have: 

For 10 7 Hz > x > 1 yr' 1 ■ Q ~ X° , with a = 0.2 - 0.4 (0.2 for partial melts) . 
For 1 yr^ 1 > x > 1CT 2 yr^ 1 : Q ~ x ■ 

For 10~ 2 yr~ x > x '■ arguably, it is still Q ~ x- ( Or maybe Q ~ l/x ?) 

Fortunately, in the study of tides in planets one never has to approach the Maxwell-time scales, 
so the controversy remaining in (j24"|) bears no relevance to our subject. In this paper we shall 
address tidal dissipation only in planets, leaving satellites for later. We shall not address the 
unique case of the Pluto-Charon resonance, nor shall we address the binary asteroids locked 
in the same resonance. Thus we shall avoid the frequency band addressed in (j23|) . We shall be 
interested solely in the frequency range described in (j22|) . It is important to emphasise that 
the positive-power scaling law (1221) is well proven not only for samples in the lab but also for 
vast seismological basins and, therefore, is universal. Hence, this law may be extended to the 
tidal friction. 

Below we provide an extremely squeezed review of the published data whence the scaling 
law ff22l was derived by the geophysicists. The list of sources will be incomplete, but a full 
picture can be restored through the further references contained in the works to be quoted 
below. For a detailed review on the topic, see Chapter 11 of the book by Karato (2007) that 
contains a systematic introduction into the theory of and experiments on attenuation in the 
mantle. 

4.2 Circumstantial evidence: attenuation in minerals. 
Laboratory measurements and some theory 

Even before the subtleties of solid-state mechanics with or without melt are brought up, the 
positive sign of the power a in the dependence Q ~ x a ma Y be anticipated on qualitative 
physical grounds. For a damped oscillator obeying 'z + 2[3z + x 2z = Q, the quality factor 
is equal to x/(2/3) , i.e., Q ~ x ■ 

Solid-state phenomena causing attenuation in the mantle may be divided into three groups: 
the point-defect mechanisms, the dislocation mechanisms, and the grain-boundary ones. 

Among the point-defect mechanisms, most important is the transient diffusional creep, i.e., 
plastic flow of vacancies, and therefore of atoms, from one grain boundary to another. The flow 
is called into being by the fact that vacancies (as well as the other point defects) have different 
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energies at grain boundaries of different orientation relative to the applied sheer stress. This 
anelasticity mechanism is wont to obey the power law Q ~ x a with a ~ 0.5 . 

Anelasticity caused by dislocation mechanisms is governed by the viscosity law Q ~ x 
valid for sufficiently low frequencies (or sufficiently high temperatures), i.e., when the viscous 
motion of dislocations is not restrained by the elastic restoring stressjj 

The grain-boundary mechanisms, too, are governed by the law Q ~ x° S though with a lower 
exponent: a ~ 0.2 — 0.3. This behaviour gradually changes to the viscous mode (a = 1) at 
higher temperatures and/or at lower frequencies, i.e., when the elastic restoring stress reduces. 

We see that in all cases the quality factor of minerals should grow with frequency. Accord- 
ingly, laboratory measurements confirm that, within the geophysically interesting band of x ? 
the quality factor behaves as Q ~ x a with a = 0.2 — 0.4 . Such measurements have been 
described in Karato & Spetzler (1990) and Karato (1998). Similar results were reported in the 
works by the team of I. Jackson - see, for example, the paper (Tan et al 1997) where numerous 
earlier publications by that group are also mentioned. 

To this we would add that in aggregates with partial melt the frequency dependence of 
Q keeps the same form, with a leaning to 0.2 - see, for example, Fontaine et al (2005) and 
references therein. 



4.3 Direct evidence: attenuation in the mantle. 
Measurements on seismological basins 

As we are interested in the attenuation of tides, we should be prepared to face the possible 
existence of mechanisms that may show themselves over very large geological structures but 
not in small samples explored in the lab. No matter whether such mechanisms exist or not, 
we would find it safer to state that the positive-power scaling law Q ~ x° j even though 
well proven in the lab, must be propped up by a direct seismological evidence gathered over 
vast zones of the mantle. Fortunately, such data are available, and for the frequency range of 
our interest these data conform well with the lab results. The low-frequency measurements, 
performed by different teams over various basins of the Earth's upper mantle, agree on the 
pivotal fact: the seismological quality factor scales as the frequency to the power of a positive 
fraction a - see, for example, Mitchell (1995), Stachnik et al (2004), Shito et al (2004), and 
further references given in these sources] 10 ! 

9 At higher frequencies or/and lower temperatures, the restoring force "pins" the defects. This leads to the 
law Q ~ (1 + t 2 x 2 )t~ 1 x~ 1 i parameter r being the relaxation time whose values considerably vary among 
different mechanisms belonging to this group. As the mantle is warm and viscous, we may ignore this caveat. 

10 So far, Figure 11 in Flanagan & Wiens (1998) is the only experimental account we know of, which 
only partially complies with the other teams' results. The figure contains two plots depicting the frequency 
dependencies of 1/ Q sheer and 1 / Q compress ■ While the behaviour of both parameters remains conventional 
down to 10 _1 Hz, the sheer attenuation surprisingly goes down when the frequency decreases to 10~ 3 Hz. 
Later, one of the Authors wrote to us that "Both P and S wave attenuation becomes greater at low frequencies. 
The trend towards lower attenuation at the lowest frequencies in Fig. 11 is not well substantiated^ (D. Wiens, 
private communication) Hence, the consensus on ([22]) stays. 
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4.4 Consequences for the tides 



4.4.1 Tidal dissipation vs seismic dissipation 

One of the basic premises of our model will be an assumption that, for terrestrial planets, the 
frequency-dependence of the Q factor of bodily tides is similar to the frequency- dependence 
( 1221 - 123]) of the seismological Q factor. This premise is based on the fact that the tidal 
attenuation in the mantle is taking place, much like the seismic attenuation, mainly due to the 
mantle's rigidity. This is a nontrivial statement because, in distinction from earthquakes, the 
damping of tides is taking place both due to rigidity and self-gravity of the planet. Modelling 
the planet with a homogeneous sphere of density p , rigidity p , surface gravity g , and radius 
R, Goldreich (1963) managed to separate the rigidity-caused and self-gravity- caused inputs 
into the overall tidal attenuation. His expression for the tidal quality factor has the form 

Q being the value that the quality factor would assume were self-gravity absent (i.e., were 
damping due to rigidity only). To get an idea of how significant the self-gravity-produced 
input could be, let us plug there the mass and radius of Mars and the rigidity of the Martian 
mantle. For the Earth mantle, p = 65 -f- 80 GPa. Judging by the absence of volcanic activity 
over the past hundred(s) of millions of years of Mars' history, the temperature of the Martian 
upper mantle is (to say the least) not higher than that of the terrestrial one. Therefore we 
may safely approximate the Martian p with the upper limit for the rigidity of the terrestrial 
mantle: p = 8 x 10 10 Pa. All in all, the relative contribution from self-gravity will look as 

— 11 3 —1 —2 23 2 

2 g p R 6 7 M 2 ^ 1 ( 6.7 x 10 m kg s ) ( 6.4 x 10 kg ) -a 
19 ~p~ ~ 70^ JR 1 * 40 (ll 10 Pa) (3.4 x 10° m) 4 * 

where 7 stands for the gravity constant. This rough estimate shows that self-gravitation 
contributes, at most, several percent into the overall count of energy losses due to tides. This 
is the reason why we extend to the tidal Q the frequency-dependence law measured for the 
seismic quality factor. 

The above estimate is conservative, because we approximated the Martian p with the 
terrestrial value of p . If the interior of Mars is cooler than that of the Earth, then its rigidity 
should be higher, and therefore the ratio (|26|) may reduce to less than one percent. 

4.4.2 Dissipation in the planet vs dissipation in the satellite 

As already agreed above, in this work we are dealing solely with the satellite-generated tides on 
the planet. The input from the planet-caused tides on the satellite will be considered elsewhere. 
Hence, here it will not be a problem for us to encounter a configuration where a satellite shows 
the same face to its primary - this kind of spin-orbit resonance will not imply a vanishing tidal 
frequency. Nor will it be a problem for us to have a satellite crossing a synchronous orbit - 
the presence of the radial term in (TlTl - [T8|) will prevent the tidal frequency from becoming nil 
at the point of crossing. The same is true for the lag related to the frequency through formula 
([3]): although the "tangential" lag 5i vanishes at such a crossing, the "radial" part of the lag 
remains finite. Hence, the overall 5 stays finite. 
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A very special situation is tidal relaxation toward the state where the planet shows the same 
side to its moon, like Pluto does to Charon0 A gradual approach to this state makes the 
tidal frequency asymptotically go to zero. Mathematically, this situation still may be tackled 
by means of (1221) until the tidal frequency x decreases to lyr -1 , and then by means of (1231) 
while x remains above the inverse Maxwell time of the planet's material. Whether the latter 
is physically satisfactory remains an open issue of a generic nature that is not related to a 
specific theory of tides or to a particular frequency dependence of Q . The generic problem is 
whether we at all may use the concept of the quality factor on approach to the Maxwell time, 
or whether we should, beginning from some low x j employ a comprehensive hydrodynamical 
model. In the current work, we shall not address this question, leaving it for another paper. 

Thus, since we are talking only about dissipation inside the planet, and are not addressing 
the exceptional Pluto-Charon case, we may safely assume the tidal frequency to always exceed 
1 yr~ x . Thence (122!) will render, for a typical satellite: 

Q ~ x a , with a = 0.2 - 0.4 . (27) 

Accordingly, (1211) will entail: 

5 ~ x~ a , with a = 0.2 - 0.4 . (28) 

4.5 The frequency and the temperature 

In the beginning of the preceding subsection we already mentioned that though the tidal 
Q differs from the seismic one, both depend upon the frequency in the same way, because 
this dependence is determined by the same physical mechanisms. This pertains also to the 
temperature dependence, which for some fundamental reason combines into one function with 
the frequency dependence. 

As explained, from the basic physical principles, by Karato (2007, 1998), the frequency 
and temperature dependencies of Q are inseparably connected. Since the quality factor is 
dimensionless, it must retain this property despite the exponential frequency dependence. This 
may be achieved only if Q is a function not of the frequency per se but of a dimensionless 
product of the frequency by the typical time of defect displacement. This time exponentially 
depends upon the activation energy A* , whence the resulting function reads as 

Q ~ [ X exp(A*/RT)] a . (29) 

For most minerals of the upper mantle, A* lies within the limits of 360 — 540 kJ mol -1 . For 
example, for dry olivine it is about 520 kJ mol -1 . 

Thus, through formulae (1291) and (121]) . the cooling rate of the planet plays a role in the 
orbital evolution of satellites: the lower the temperature, the higher the quality factor and, 
thereby, the smaller the lag 5 . For the sake of a crude estimate, assume that most of the tidal 
attenuation is taking place in some layer, for which an average temperature T and an average 
activation energy A* may be introduced. Then from (T2"9"j) we have: AQ/Q w a A* / RT . For 
a reasonable choice of values a = 0.3 and A* = 5.4 x 10 5 J/mol , a drop of the temperature 
from T = 2000 K down by AT = 200 if will result in AQ/Q « 1 . So a 10% decrease of 
the temperature can result in an about 100 % growth of the quality factor. 

Below we shall concentrate on the frequency dependence solely. 

11 Such a complete locking is typical also for binary asteroids. Since at present most asteroids are presumed 
loosely connected, and since we do not expect the dependencies (|22l - [23|) to hold for such aggregates, our theory 
should not, without some alterations, be applied to binaries. 
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5 Formulae 



The tidal potential perturbation at the surface of a planet due to a point satellite is given by 



W{i) = J2 ^(7. 



(30) 



n=2 



where 



so that 



W 



Gm R 2 



GmR n 
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P n (cosj) 



( 3 cos 2 7 — l)+ — (5 cos 3 7—3 cos 7 ) + 



(31) 



(32) 



R being the mean equatorial radius of the primary, and the angle 7 being reckoned from 
the planetocentric vector pointing at the satellite. Being small, the tidal perturbations of the 
potential yield a linear response of the planet's shape, deformations proportional to the Love 
numbers. These deformations, in their turn, amend the external potential of the planet with 
an addition, which would look as 



U{l) 



x>(*) 

n=2 v ' 



n+1 



W n ( 7 , r) 



(33) 



had the planet's response to the satellite's pull been synchronous. In reality the bulge lags, as 
if it emerges beneath a fictitious satellite lagging relative to the real one. If the real satellite is 
located at r = r(u), where v is the true anomaly, then the fictitious satellite is at 



r f = r+f , / = At (v p x r - v) 
In these terms, the external potential of the planet will read: 



(34) 



n+l 



U{Si 



n=2 \ f 



W n (-y, r) 



(35) 



the angle b~\ = arccos -> 4 being reckoned from the bulge. (Mind that the lag 5 is, generally, 

|y| [T"jf I 

different from 8\ , and coincides with it only in the limit of vanishing e . ) To calculate the 
back-reaction exerted by the bulge on its progenitor satellite, one should equate 5\ and 7 . 
Then (1331) will assume the form: 
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= ife ( 3 (r r ?) 2 - r]r 2 ) + (5 (f> • r) 2 - 3r 2 r 2 ) + ... (36) 

where r = \r\ and rf = , while the constants are given by 

k 2 GmR 5 . k 3 GmR 7 
^2 = g ' As = 2 ' " ' ' ( } 

To (131)]) we must add the potential due to the tidal deformation of the satellite by the planet. 
That input contributes primarily to the radial component of the tidal force exerted on the 
moon, and redounds to the decrease of eccentricity (MacDonald 1964). Here we omit this 
term, for our goal is to clarify the frequency dependence of the lag. In a subsequent work, 
devoted to orbit integration, we shall take that term into account. 

Another important comment worth mentioning is that in practical applications one often 
has to take into account the k$ and sometimes even the k& terms of the expansion for the tidal 
potential. (Bills et al 2005) Once again, since in the current paper our goal is to explain the 
new model and not to perform whatever numerical integration, then we shall restrict ourselves, 
for brevity, to the leading term proportional to k 2 ■ 

— » 

Our intention is take into consideration the frequency-dependence of the lag J 7 , but not of 
the parameter k 2 . While the dependence J~(x) will be derived through the interconnection 

— » 

of J- with 6 and therefore with Q(x) > the value of k 2 will be kept constant. This can be 
justified by means of the following formula derived by Darwin (1908) under for a Maxwell body 
(see also Correia & Laskar 2003): 



Mx) = k fluid 



i + x 2 rflu 2 



1 + (x 2 V 2 /^ 2 ) (1 + l9fi/(2gpR)) 2 



kfiuid being the fluid Love number. This is the value that k 2 would have assumed had the 
planet consisted of a perfect fluid with the same mass distribution as the Maxwell-body planet 
considered. Notations \i , p , g , and g stand for the rigidity, mean density, surface gravity, 
and the radius of the planet. The letter rj signifies the viscosity. For a terrestrial planet's 
mantle, the latter will be close to 10 22 kg/(m- s), give or take an order or two of magnitude. 
For the other parameters we may keep using the estimates from subsection 4.4.1. This will 
entail: x 2 V 2 /^ 2 = (x ' 10 11 s) 2 , whence we see that in all realistic situations pertaining 
to terrestrial planets the frequency-dependence in Darwin's formula will cancel out. Thus we 
shall neglect the frequency-dependence of the Love number k 2 (but shall at the same time 
take into account the frequency- dependence of Q , for it will induce frequency-dependence of 
all three lags). 

This said, we can apply the standard rule of thumb according to which calculation of the 
tidal force consists of the following four steps: 

1. Differentiate fl36|) with respect to r , keeping r) constant. 
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2. After the differentiation^ 2 ] insert the expression rf = r+ f . 

— * 

3. Substitute / with the product At ( uj p x r — v) . 

4. Express everything in terms of the orbital elements. 



This will entail an unaveraged expression for the tidal force in terms of the satellite's or- 
bital elements. It can then be plugged into the Gauss planetary equations. Averaging thereof 
will yield the equations describing the evolution of the secular parts of the orbital variables. 
Steps 1 and 2 lead us to the following formula for the force: 
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(38) 



taken that 0(f'~/r 2 ) = 0(5 2 ) < 1 , and assuming that the satellite is pulled by the bulge 
it itself has created, with no other moons involved. 

To perform step 3, one should use connection ffT9~]) between the position and angular lags: 



or, equivalently: 



where 
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(41) 



is the unit vector pointing in the lag direction, and 5 is a known function of the frequency x ■ 
Be mindful that all equations (1381 - 141!) are exact, in that they do not make use of our (or 
any other) specific model of tides. It is of a crucial importance that the subtended angle (the 
"tangential lag") 5i does not appear in (!38l) . The expression for the force depends solely upon 

— * 

f , which in its turn is interconnected with the total lag S and the quantity x{ v ) ~ simply 
because 5 is defined as 5 = \f\/r = r~ l \uj x r — v\ At and because x{ v ) is defined as 
2r~ l \uj x r — v\ . We wish to emphasise that in the end of the day the actual physical force 

— » 

J- depends not upon the "tangential lag" b\ but upon the total lag 5 . 

The realm of our particular model begins when, for the reasons explained in subsection 2.4, 
we interpret xiy) as the instantaneous tidal frequency, in the spirit of the WKB approximation. 
This immediately brings into play the instantaneous quality factor Q as a function of the 



12 At this point wc are on the horns of dilemma. Should we first take the partial derivative an( i 

impose the constraint r/ = r + / afterwards, or should we perform these operations in the opposite order? 
The former option is conventional and, fortunately, correct, though its justification requires some effort. 
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instantaneous frequency: Q = Q{x{ v )) ■ Our knowledge of the scaling law (J271 - will then 
enable us to proceed within the framework of our model. Since, according to (1281) . 5 ~ x~ a > 



then (I2T]) necessitates for the quality factor: Q 
conveniently expressed as 

Q = £ a X a , 



This power scaling law may be 



£ a simply being the dimensional factor emerging in the relation Q ~ . As mentioned in 
subsection 4.5, cooling of the planet should become a part of long-term orbital integration. It 
enters the integration via evolution of this factor £ a . If we suppose that most of the tidal 
attenuation is happening in some particular layer, for which an average activation energy A* 
and an average temperature T may be defined, then (1291) will result in: 
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T Q being the temperature of this layer at some fiducial epoch. Then the physical meaning 
of the parameter £ becomes evident: if the planet were made of a homogeneous medium, 
with a uniform temperature distribution, and if dissipation in this medium were caused by one 
particular physical mechanism, then £ would be a relaxation time scale appropriate to this 
mechanism (say, the time of defect displacement). For an actual celestial body, £ may be 
interpreted as a relaxation time averaged (in the sense of Q = £ a x a ) over the body's layers 
and over the various attenuation mechanisms acting within these layers. This way, £ is an 
integral parameter describing the overall tidal dissipation in an inhomogeneous body. 
As 022]) yields 5 « 1/(2 Q) = (l/2)£- a X ~ a , then ([33-1101) S iv e for the position lag: 

/ = r 6 f = - r £~ a X ~ a f , (42) 



and for the time lag: 
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£ being the afore introduced integral parameter, and x being a known function (1181) of the 
orbital variables. All in all, we have: 
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Formulae (|46p . (|44|) . and (|38p are sufficient to both compute the orbit evolution and trace the 
changes of the time lag. 

In neglect of the 0(e 2 ) terms, (|4"4"1 - |4U|) simplify to 

-j /„ N~ a a uj v x r — v , _ N 

f " < £ *> TT7^ \lxr--vi • (47 > 



X~2y [n(l + 2e cosz/) — u^cosz] + cj^ sin 2 z cos 2 (u; + v) , (48) 

and 

-(a+l)/2 

-(a+l) -a/ 2 2 2 2 \ 

At ^ 2 £ I [ra (1 + 2e cosi/) — aj p cosz] + c<j p sin i cos (uj + i/) ) . (49) 

In the MacDonald-Kaula-Gerstenkorn theory the geometric lag angle 5 is constant (i.e., a = 
0). The Singer-Mignard theory suggests that we treat the time lag At as a constant (i.e., 
set a = —1). As explained in the preceding sections, neither of these choices corresponds 
to the geophysical realities. The physics of attenuation dictates to us that for a realistic 
terrestrial planet, the exponential a lies within the range of 0.2 — 0.4 , except in the case of 
near-synchronous satellites for which it increases to unity. 

To write down the planetary equations in the Euler-Gauss form, either the exact expression 
( 144|) or its approximation (j4Tj) should be plugged into (l38]) f^l The equations can then be 
integrated. One possibility would be a straightforward numerical integration. An alternative 
option would be to perform analytical averaging over the true anomaly v , and only afterwards 
to integrate the averaged Euler-Gauss-type equations describing the evolution of the secular 
parts of the satellite's orbital elements. 

A numerical simulation of the Mars-Phobos dynamics, in neglect of the primary's cooling 
and in neglect of the tides on the satellite, will be presented in our next publication (Efroimsky 
& Lainey 2007b). Some preliminary results of that study have been reported in our talk 
Efroimsky & Lainey (2007a). According to those computations, Phobos will fall on Mars in 45 
Myr from now. This estimate is 50 percent longer than the one following from the tidal models 
used in the past. An overwhelming share of this difference comes from the actual frequency 
dependence of attenuation, and only a small part comes from our usage of the instantaneous 
tidal frequency instead of the Kaula expansion. This demonstrates that the currently accepted 
time scales of dynamical evolution should be reexamined using the actual frequency dependence 
of the lags. 

13 Equivalently, one can simply plug ((2]) into (|38|) . to get 

T = - 3 fc2 G ™ R [ff r 2 _|_ (2f* (f-v) + r 2 (fx (2 + if)) At] + ... , 

and can insert it into the Euler-Gauss-type planetary equations. Formally, this will render the same unaveraged 
equations for the elements' evolution, as the equations derived by Mignard (1979, 1980). However, integration of 
our equations should be carried out with an important proviso that now At is not a constant but is expressed 
by (TJ3]), i.e., depends (through (TJS]) or (T49")0 upon the current values of the orbital elements and the true 
anomaly. This will cardinally alter the results of integration. 
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6 Conclusions 



In our model, we operate with the overall quality factor Q that describes both components of 
the tidal friction - those associated with the "tangential" and "radial" lags. By the "tangential" 
lag we understand the angle 8± subtended at the planet's centre between the satellite and the 
bulge. It emerges because the location of the bulge lags. The "radial" lag is called into being by 
the fact that the height of the bulge lags too. (This latter lag exists only for noncircular orbits, 
and vanishes after an orbit gets circularised.) Our model is based on the observation that the 
overall quality factor Q is interconnected, via the standard relation (120|) . not with the "tan- 



gential lag" 81 , but with the total lag 8 introduced as y (tangential lag) + (radial lag ) 



The model also employs an instantaneous tidal frequency x > which is a function of the true 
anomaly v of the satellite, - a concept vaguely prompted by the WKB method in quantum 
mechanics. We define this instantaneous tidal frequency x as the quantity, insertion whereof 
in ( jTTi) or ( |T2l) gives the correct dissipation rate. This definition yields the expression ( |T8l) 
for x through the satellite's orbital elements and the true anomaly. For a circular equatorial 
orbit, the right-hand side of that formula reduces to the well known expression 2\n — uj p \ . 

So defined total lag 8 , instantaneous frequency x > an d the overall quality factor Q are 
always (i.e., for an arbitrary eccentricity and inclination of the moon) connected among them- 
selves by the same relations - |2"01 as in the trivial case of an equatorial circular orbit. 
This has enabled us to circumvent the difficulties of the standard, Kaula-expansion-based tidal 
theory by completely avoiding whatever decomposition of the tide. The instantaneous tidal 
frequency x{ v ) an d the instantaneous overall quality factor Q(x) — Q( x{ v ) ) are t ne on ly 
frequency and the only Q factor entering the solution. This model is the only means of escape 
in nonlinear situations where the quality factors of the harmonics are badly defined (i.e., when 
the "Goldreich's admonition" is in force - see subsection 2.2 above). 

Within our model we permit the overall quality factor and, therefore, both the overall 
angular lag S and the time lag At to depend upon the instantaneous tidal frequency x ■ While 
in the Gerstenkorn-MacDonald-Kaula theory the geometric lag is postulated to be constant, 
in the Singer-Mignard theory it is the time lag that is assumed constant. However, neither of 
these two choices conform to the geophysical data within the frequency band of our interest. As 
the overall angular lag S is inversely proportional to the tidal Q factor, the actual frequency- 
dependence of both 5 and At is unambiguously defined by the frequency-dependence of Q . 

The quality factor is wont, according to numerous studies, to obey the law Q ~ x a > where 
a lies within 0.2 — 0.4 . For this reason, the time lag At , employed in the Singer-Mignard 
theory as a fixed parameter, in reality is not a constant but a function (??) of the frequency. 
Hence, through (|4"6T) . At becomes a function of the current values of the orbital elements of 
the satellite, and so becomes 8. This clause alters the integration of the planetary equations. 

In a subsequent publication we shall deal with the numerics, in application to the possible 
history of the Martian satellites (Efroimsky & Lainey 2007b). 
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Appendix 



A.l. The expression for i | Q v x r — v\ 

Our goal is to prove formula (TTTT) . To that end, let us start out with the evident equality 



1 , - 



r 2 I ~P 



u)„ x r - v\ 2 = \ [ (ti p x rf + 2 - 2 (<3 p x f) ■ v] . (50) 



The satellite velocity v can be expanded into its radial and tangential components: v r — f 
and v = r v , with v standing for the true anomaly. Thence we get: 



' 1 - "* " |2 - \ [(£ p xr) 2 + r 2 u 2 + v 2 - 2 (cD p xr)-^] (51) 



— \u> p x r — v 



dM dv d f a (1 — e ) \ (1 + e cos v) a (1 — e ) e sin n a e 



where 



v = 7* = — i i = ti = — — sin v 

dt dM dv \l + e cosv J (l-e 2 ) 3/2 (1 + ecosz/) 2 y/l - e 2 

To continue, we need expressions for ( uj p x f) 2 and — 2 ( u5 p x r ) • v . To derive these, we 
introduce a nonrotating planetocentric Cartesian coordinate system (x, y, z) with z pointing 
northward along the planet spin axis, and with x , y lying within the equatorial plane. (In our 
model, we neglect the planetary precession, nutation, Chandler wobble, or the polar wander, 
so this system may be referred to as quasi- inert ial.) The satellite's longitude of the node may 
be reckoned, say, from the axis x .) In this coordinate system, the planet's angular spin rate is 



while the satellite position 



has components: 



uj p = uj p z , (53) 



xx + yy + zz (54) 



1 - e 2 

x = a f cos f2 cos (u + v) — sin Vt sin (u + u) cos i } , (55) 

1 + e cosu 1 1 

1 - e 2 

y = a [ sin cos (u + u) + cos sin (u + v) cos i ] , (56) 

1 + e cosu 



1 - e 2 

z = a sin (u + v) sin % . (57) 

1 + e cos v 



Formulae ([53] - [54]) entail: 

w p x f = — Up y ~k + Up x y , (58) 
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whence we get the desired expression for (uj p x r ) 2 : 

(£ p x r) 2 = uo 2 p (x 2 + y 2 ) = u 2 p (r 2 - z 2 ) = u 2 p r 2 ( 1 - sin 2 i sin 2 (cj + u) ) . (59) 
The satellite velocity 

v = v x x + Vy y + v x z (60) 

has Cartesian components 

[ — cosf2 sin(u; + u) — sinf2 cos (a; + v) cos i + 



e ( — cosfi sincj — sinfi coscj cosi) 



[ — sinf2 sin(u; + z/) + cosfi cos(c<j + i/) cos % + 



e ( — sin f2 sin + cos fl cos a? cos z 



Together, ( 1581 and (ED - [63]) will yield the needed expression for — 2 ( uj p x r ) ■ v 



and 



(61) 



(62) 



= sin i [ cos(cj + v) + e coscj] , (63) 



— 2 (uj p x r) ■ v = — 2 uj p a 2 n yl — e 2 cos i (64) 
Finally, recall that the planetocentric distance r to the satellite is 

1 - e 2 

r = a — , 65 

1 + e cos v 

wherefrom the following auxiliary formulae ensue: 

-2{uj p xr)-v (1 + ecosz/) 2 . . . 

^ = -2fiw p jj^— cos z = — 2 uo p v cos z (66) 

r (1 — e 2 ) ' 



f 2 / n a e . \ 2 ( a (1 — e 2 ) \ 2 22/-, 2\~ 3 • 2 /-, \ 2 /^\ 

-4 = ; sin z/ — i - = n 2 e 2 (1 - e 2 ) sm 2 v (1 + e cos z/) . (67) 

Plugging of (J52D, dHSD, and ([61 -[65]) into ([SI]) then leads to: 
\ (jj v x r - v = q 



[u 2 p r 2 ( 1 - sin 2 i sin 2 (cj + z/) ) + r 2 z> 2 + w 2 - 2 (u3 p xr) -v] 



oo 2 p cos 2 z + oo 2 p sin 2 z cos 2 (cj + z/) + z> 2 + (— ) 2 + 2 ^ ^ r ^ - , (68) 
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while the subsequent insertion of ([HS] - EZD will bring us to 

— | uj p x r — v\ 2 = ujp cos 2 i + Up sin 2 i cos 2 (u + v) + z> 2 + — 2u) p v cos i 



= ( cj p cos z — z> ) + cj p sin i cos (c<j + v) + ^ — J . (69) 

This proves formula ffl7|) . 

One can introduce the mean principal frequency by averaging % over one revolution of the 
satellite. Applying the averaging rule (271)- 1 j o 27r dM = ^vr)" 1 (1-e 2 ) 3 / 2 Jf(l+ e cosu)~ 2 du 
to the above expression for \ > we arrive at 



/ 1 5 

(x) = 2 W (n — ujp cosz) 2 + - ujp sin 2 z + - ri 2 e 2 + 0(sin 4 z) + 0(e 2 sin 2 i) + 0(e 4 ) 



A. 2. The quality factor and the angular lag. 

The goal of this section is threefold. First, we wish to remind to the reader why, in the 
first approximation, the quality factor Q is inversely proportional to the phase lag. Second, 
we shall remind why the phase lag is twice the geometric lag angle, as in formula (12"U1) above. 
Once again, while a comprehensive mathematical derivation of this fact can be found elsewhere 
(see the unnumbered formula between equations (29) and (30) on page 673 in Kaula (1964)), 
our goal here is to illustrate this, very counterintuitive result, by using the simplest setting. 
Third, we shall explain why in general (i.e., for nonequatorial noncircular orbits) the Q factor 
is interconnected not just with the "tangential lag" (the subtended geometrical angle 5\ on 

I 2 2~ 

Fig. 1) but with the total lag 5 defined as y (tangential lag) + (radial lag) . (While the 
"tangential" lag 5\ emerges because the position of the bulge lags, the "radial" lag comes into 
play due to the fact that also the height of the bulge lags.) 



A. 2.1 The case of a circular equatorial orbit. 

We shall begin with the simple case of an equatorial moon on a circular orbit. At each point 
of the planet, the tidal potential produced by this moon will read 

W = W cosxt , (71) 

the instantaneous frequency being given by 

X = 2 \n - u p \ . (72) 

Let g denote the free-fall acceleration. An element of the planet's volume lying beneath the 
satellite's trajectory will then experience a vertical elevation of 

C = — cos( X t - 25,) . (73) 
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Accordingly, the vertical velocity of this element of the planet's volume will amount to 

u = C — — X — ° srn (x^ ~~ 2 5) = — x — - (sinx^ cos 25, — cosxt sin 2(5,) . (74) 

g ' g 

The expression for the velocity has such a simple form because in this case the instantaneous 
frequency % ls constant The satellite generates two bulges - on the facing and opposite 
sides of the planet - so each point of the surface is uplifted twice through a cycle. This entails 
the factor of two in the expressions (1721 for the frequency. The phase in (1751) . too, is doubled, 
though the necessity of this is less evidential 

The energy dissipated over a time cycle T = 2ir/(x) , per unit mass, will, in neglect of 
horizontal displacements^ be 

f T ( 9W\ , ( W a \ dW D f t=T nr . n . 

&E cycle = J it I - — \dt= - I -X— 1 -q^t J cosxt {smxt coa26 1 - cos^t sin 2^ 



— x k — sin 25, — / cos \t d(yt) = 7: — tt sin 25 

g or x J x t = o g dr 



14 Equally simple would be a hypothetical case of an inclined circular orbit about a slowly spinning primary 
(tup <ti n). Here the frequency (|18p would, for an arbitrary inclination, be approximated simply with n, i.e, 
with a constant. Hence all derivations presented in this subsection would remain in force. 

15 Let x signify a position along the equatorial circumference of the planet. In the absence of lag, the radial 
elevation at a point x would be: 

( = — - cosk(x — v t) , v = R<j , 
g 

v being the velocity of the satellite's projection on the ground, R being the planet's radius, and a being simply 
|n — uj p \ because we are dealing with a circular equatorial orbit. The value of k must satisfy 

k v = 2 a , i.e., k v = x > 

to make sure that at each x the ground elevates twice per an orbital cycle. The above two formulae yield: 

k R = 2 . 

In the presence of lag, all above stays in force, except that the formula for radial elevation will read: 

W 

( = cos k(x — v t + D) , where D = R 8 1 , 

D being the linear lag, and 5 1 being the angular one. Since kv — 2 , we get: 

cos [k (x — v t + R S 1 ) ] = cos [k x — k v t + k R 5} — cos [k x — (k v t — 2 S ± )] , 
so that, at some fixed point (say, at x = ) the elevation becomes: 

((t) = ^ cos(fc vt - 28.) . 
g 

We see that, while the geometric lag is S 1 , the phase lag is double thereof. 

16 In Appendix A. 2. 3 below, it will be explained why the horizontal displacements may be neglected. 
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while the peak energy stored in the system during the cycle will read: 



E 



pea k 



T/4 ( dW\ , / W \ dW f t=T/4 nr . nr ,_, 

u I — ^— I at = — I — % — I — — / cosx^ (sin cos2o 1 — cosx^ sin 2^ J at 
n \ or J \ g J Or J t=0 



W dW 

2(7 x~ 

g or 



cos 25. /-^ = ^ 2 r/ 

cosxt smx* a(x*J - 



X J x t = o 



811125, f xt = 7T/2 2 . 
/ cos a(xt) 

X <V = o 



g <9r 



cos 25, 



7T 



sin 2^ 



(76) 



whence 



Q- 1 



AE 



cycle 



2ttE 



pea k 



7r sin 2(5, 



27r ^ COS25, 



| sin 2^ 



tan 25, 



(77) 



A. 2. 2 The general case of a nonzero eccentricity and inclination. 

In this subsection, we shall switch from the specific case of an equatorial circular orbit to 
the generic case of an inclined elongated one. If we want to keep using the ansatz (!72|) . we 
must now accept that the tidal frequency x is n °t a constant but a variable quantity x( v ) 
or, if we choose to parameterise the motion with time instead of the true anomaly, x{t) ■ 

Now we have the following expressions for the tidal potential and the (delayed by At ) 
vertical elevation: 



W = W Q cos[x(*) t 



(78) 



C = — cos{ X (t) [t - At}} 



(79) 



Back in (1T3]) . the quantity x At had an evident meaning of the doubled "tangential lag" angle, 
so we preferred to write it as 25 1 . Here, though, x(t) At lacks such a straightforward meaning, 
and therefore we would prefer to deal with At . 

Generalising (1751) . we must keep in mind that \ is no longer constant: 



AE 



cycle 



*x(*) = 2w dw \ Wo ( dW \ /^M = 2 - 

tx(t) = o \ dr J g V dr J J tx(t) = 



cos[x(t) t] dcos\x(t)t - x(t)At_ 



W dW f tx{t) = 2n r , w . d cos [ X (t)t - X (t)At] d[ X (t)t - X (t)At] 

^ — / COS Y t t T7 — j-r . , . , — — r — -r- s ; fit Ylt) t 

g dr J tx(t)=0 d[ x {t)t - x {t)At] d[ X (t)t] 



W dW /"*W = 2w 
g dr J tx(t) = 



d[x{t)At] dt 



dt 



d[ X (t)t) 



cos [x(t)t] sin [ X {t)t- X (t) At ] d [ X (t) t } 
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Wo dWo 

p. Q r 



rtx(t) = 2w 

/ [1 - 0(xAt/x)} (cosxt cos x At - sinxt sin x At) sinxtd(x^) • ( 

Jt x (t) = o 



Let us now recall that x{t) is a slow function of time, i.e., that x/x 2 can be treated as a 
small dimensionless parameter Also keep in mind that At, too, is small (i.e., that x^t 
too can be treated as a small dimensionless parameter). It will then be legitimate to neglect, 
in the square brackets, the term 0(xAt/x) = O ( (x At) (x/x 2 ) ) whose order of smallness is 
the same as that of the product of these two small parameters. For the same reason, in the 
right-hand side of flHUj) one can approximate cosxAt with unity and can treat sinxAt as a 
constant. Thence we shall obtain an analogue to ( [75]) : 

W dWo 

dr 

In the same approximation of x(t) being a slow function, we can get an analogue to ( 1761) : 



= " -i^ * sin [x(t) At] . (81) 



AE „ = 

peak Q 



i cos [x(t) At] - | sin [x(t) At] 



52) 



Both in the latter and in the former formulae integration over the entire volume of the body 
is implied@ From (1HT1) and ( 1821) . the following generalisation of ( 1771) ensues: 

Q" 1 = " = -1 n ^sin[x(t)At] tan[x(0A£] . (83 ) 



peak 



% cos [x(t) At] - | sin [x(t) At] 



So in the general case the factor Q is connected not to the subtended angle but to x{t) At/2 . 
A. 2. 3 The final touch 

In our above derivation of the interrelation between Q and 5 , we greatly simplified the 
situation, taking into account only the vertical displacement of the planetary surface, in re- 
sponse to the moon's pull. Here we shall demonstrate that this approximation is legitimate, at 
least in the case when the planet is modelled with an incompressible (though not necessarily 
homogeneous) medium. Thence we shall see that the calculations offered in a simplified form 
in A. 2.1 and A. 2. 2 can be expressed in the fancy language of continuum mechanics. 

As a starting point, recall that the tidal attenuation rate within the planet is well approx- 
imated with the work performed on it by the satellite 

The attenuation rate within a tidally distorted planet is: 

E = - J P V VW d 3 x (84) 

— * 

where p , V and W are the density, velocity, and tidal potential inside the planet. To simplify 
this expression, we shall, following Goldreich (1963), employ the equality 

p v VW = V • (p V W) - W V ■ Vp - W V • (p V ) = V • (p V VW) - W V ■ Vp + W -k- 

ot 



It is straightforward from (Ti8]) that x/x 2 = O ( e{n/uj p ) 2 ) + O ( {n/u p ) sin 2 i ) . 

18 Just as in the previous section, here we have in fact performed all calculations for an elementary volume, 
and did not bother, in (|8Tj) or ([82]) . to integrate over the entire volume. We afford this omission, because this 
integration would not affect the 5±— dependent multipliers and therefore would yield the same functional form 
for the dependence Q(x(t) A) ■ The reason for this is explained in Appendix A. 3. 

19 A small share of this work is continuously being spent for decelerating the Earth rotation. 
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For a homogeneous and incompressible primary, both the V WVp and dp/dt terms are nil, 
whence 

E = - J pW V ■ nd 3 x , (86) 

n being the outward normal to the surface of the planet. We immediately see that, within the 
hydrodynamical model, it is only the radial elevation rate that matters. 

Now write the potential as W = W Q cos(x(i) t) , where x(t) is a s l° w ( m the sense ex- 
plained in A. 2. 3) function of time. Since the response is delayed by At , the surface-inequality 
rate will evolve as Vn ~ sin{x(t) [t — At] } . All the rest is as in subsection A. 2. 2 above. 



A. 3. Our calculations in the light of Zschau's theory 

Amendments to the conventional theory of bodily tides were offered by Zschau (1978) whose 
starting point was that attenuation of land tides, in distinction from earthquakes, takes place 
both due to rigidity and self-gravity of the body. 

For this reason, Zschau introduced not one but two different lags. The first was the ordinary 
lag 5 depicting the delay between the bulge on the planet's surface and the position of the 
tide-raising satellite. (Zschau called this quantity an •) The second was the (neglected in the 
literature thitherto) lag of the secondary potential emerging due to the planet's deformation. 
(Zschau called it ax •) These two lags are, generally, different from one another, and this 
difference reflects the fact that the phase delay between the tidal stress and strain varies over 
the depth of the planet and is, generally, different from that on the planet's surface. 

The traditional approach (in whose spirit our Appendix 2 above is written) implies that, as 
long the relation Q^ 1 = tan(25) is proven for an elementary volume, it can be trivially applied 
also to the entire volume of the body (provided the body is assumed homogeneous). In the 
theory of Zschau, this generalisation is no longer trivial, even if the body is homogeneous. This 
circumstance entails, for the overall quality factor of the planet, an expression more complicated 
than Q^ 1 = tan 25 . According to Zschau (1978), the inverse quality factor should read as 

Q' 1 = - \ [(1 + k) tan2a H - k tan2a x ] , (87) 

k and h being some Love numbers. Zschau emphasises, though, that the second term in the 
square brackets becomes irrelevant when k <C 1, i.e., when the self-gravity may be neglected. 

Fortunately, it turns out that the self-gravity does not bring much into the overall dissi- 
pation process, because the mantle damps tidal displacements, much like the seismic waves, 
mainly due to its rigidity. If we model the body with a homogeneous sphere of density p , 
rigidity p , surface gravity g, and radius R, then the estimate presented in subsection 4.4.1 
above tells us that self-gravitation contributes, at most, several percent into the overall energy 
dissipation due to tides. 

This fact entails two important consequences. One is that we may extend to the tidal Q 
the frequency-dependence law measured for the seismic quality factor. The other consequence 
is that the extra term —k tan2a^ in Zschau's expression for Q' 1 may be ignored (unless we 
are talking about some hypothetical terrestrial exoplanet several times larger than Mars). As 
soon as this term is omitted, 

Q' 1 = - 1 , tan 25 
1 — h 
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where we are using the conventional notation 5 instead of Zschau's an . Effective renormali- 
sation of Q by a factor of 1/(1 — h) is taking place due to gravitational prestressing, and is 
of a marginal importance to our study. We see that for the terrestrial planets (and, the more 
so, for small bodies) one still may use the conventional simple relation between the energy 
dissipated in one cycle divided by the peak energy and the tidal lag. 
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